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P-] ■ Understanding the causes and effects of spatial aggregation is one of the most 

, fundamental problems in ecology. Aggregation is an emergent phenomenon arising 

from the interactions between the individuals of the population, able to sense only 
-at most- local densities of their cohorts. Thus, taking into account the individual- 
q-( level interactions and fluctuations is essential to reach a correct description of the 

population. Classic deterministic equations are suitable to describe some aspects 
of the population, but leave out features related to the stochasticity inherent to the 
discreteness of the individuals. Stochastic equations for the population do account 
for these fluctuation-generated effects by means of demographic noise terms but, 
^T) ' owing to their complexity, they can be difficult (or, at times, impossible) to deal 

with. Even when they can be written in a simple form, they are still difficult to 
I/"") \ numerically integrate due to the presence of the "square-root" intrinsic noise. In 

this paper, we discuss a simple way to add the effect of demographic stochastic- 
ity to three classic, deterministic ecological examples where aggregation plays an 
important role. We study the resulting equations using a recently-introduced inte- 
gration scheme especially devised to integrate numerically stochastic equations with 
demographic noise. Aimed at scrutinizing the ability of these stochastic examples 
to show aggregation, we find that the three systems not only show patchy configu- 
rations, but also undergo a phase transition belonging to the directed percolation 
universality class. 

1 Introduction 

The connection between different scales of observation is a central issue in many dis- 
ciplines. The most suitable level of resolution to scrutinize a given system depends 
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essentially on the type of questions to be answered. A description at the level of the indi- 
vidual components -in which all elementary interactions are taken into account- contains 
all the information, but at the price of being, in general, intractable analytically or, in 
some cases, even computationally. In the same way as in physics the Navier-Stokes equa- 
tions are used to describe fluid dynamics in terms of coarse-grained fields (by-passing the 
description in terms of individual molecules), in ecology it is usual to resort to population- 
level dynamical equations encapsulating the most relevant features of groups composed 
of a large number of individuals. In this approach, individuals are replaced by "fields", 
which account for the density of organisms at specific points of space and time. Thus, 
deducing equations for those density fields starting from individual interactions, without 
losing relevant information, is an essential task. 

Almost two decades ago, Durrett and Levin reviewed the standard (coarse-grained) 
modeling approaches in population ecology pp. They discussed different levels of de- 
scription (individual level and continuous equations, in both their spatially explicit and 
implicit versions) as well as their mutual interconnections. Each of these approaches is 
able to reproduce some features correctly, even if with limitations. All the continuous 
descriptions presented in [1] are deterministic, and therefore suppress the stochasticity 
inherent in the individual level. Even the "hydrodynamic approach", introduced in [1] 
and extended in [2], which does account for the discreteness of the individuals, eventually 
leads to deterministic continuous equations at the population level, in which stochasticity 
is averaged away by assuming a Poissonian distribution for the number of individuals. 

In most cases, deterministic descriptions suffice to capture successfully the population- 
level phenomenology and are likely to be analytically tractable. Thus, they may provide 
a global understanding of the system phenomenology and allow us to quantify the effect 
of different factors. However, in some other cases neglecting stochastic effects leads to the 
loss of relevant information. In particular, deterministic approaches fail to describe noise- 
induced effects, which can be crucial in certain situations especially for low population 
densities and low spatial dimensions |3j. To focus our presentation, here we restrict 
ourselves to the role of noise in the problem of pattern formation in ecology. We shall 
discuss three different examples. 

Based on the early work of Turing [4] and others [5], Levin and Segel (and, indepen- 
dently, Okubo [S]) introduced a model describing the dynamics of interdependent phy- 
toplankton and zooplankton populations at a deterministic (highly coarse-grained) level 
[7j. The corresponding differential equations develop characteristic (Turing) patterns of 
aggregation when the differential diffusion of the two species is large, while homogeneous 
stationary states are obtained otherwise [7]. In particular, the region in the parameter 
space for which patterns are observed is relatively narrow as opposed to real planktonic 
populations, which typically appear in patchy distributions even at scales not traceable 
to physical forcing. Moreover, the assumption of passive diffusion for zooplankton cannot 
be biologically justified, since zooplankton can actively aggregate, and this is reflected in 
the fact that model predictions of phytoplankton being more patchily distributed than 
zooplankton are not borne out in empirical data. 

Aimed at clarifying this paradoxical situation, Steele and Henderson considered a 
model very similar to the Levin-Segel Model (LSM), but including stochasticity (i.e. a 
Gaussian white noise) in a parameter value [8]. Such a noise term -mimicking intraspe- 
cific variability in either phytoplankton or zooplankton- allows one to scrutinize whether 
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ecological interactions suffice to explain patchiness, with no need to invoke external (envi- 
ronmental) variation or fluctuations, or active aggregation. Actually, the introduction of 
noise results in a widening of the region of patchiness in the parameter space and, there- 
fore, serves as a possible explanation for the origin of pattern formation in some real-life 
examples [8]. In the same spirit, Butler and Goldenfeld recently showed that, indeed, 
fluctuations can expand the pattern region in the LSM [9]. First, they showed that the 
naive method of adding by hand a simple white noise to the deterministic LSM is able to 
induce patterns of aggregation under less constrained condition. Then, they introduced 
an individual- level model representing the interactions between the two planktonic types, 
for which they performed a rigorous scaling-up by employing standard techniques from 
statistical physics and field theory [HJ QUI EH]- This approach allowed them to derive a 
set of stochastic (Langevin) equations whose deterministic part coincides with that of the 
original LSM but where, additionally, there is a non-trivial stochastic part or noise (that 
includes non-trivial correlations, cross-correlations, and diffusive noise) [9], [10]. These de- 
mographic noise terms are a direct consequence of the stochastic nature of the underlying 
birth and death processes. Their elegant analytical calculation reveals that demographic 
or "intrinsic" noise greatly enlarges the region of the parameter space where pattern 
formation occurs. Hence, even in the absence of either environmental or intraspecific 
variability, patterns appear much more generically than in the purely deterministic LSM 
as a mere consequence of demographic noise. 

A second example in which intrinsic stochasticity due to discreteness plays a relevant 
role in pattern formation is the study of vegetation growth in semi-arid environments. 
Klausmeier introduced a continuous model consisting of two coupled differential equations 
representing the interactions between vegetation biomass and water [12] . In the presence 
of terrain slopes or soil inhomogeneities, the model develops Turing and/or disordered 
patterns [121 H3] . However, the model is not able to generate patterns in plain terrains, 
for which only homogeneous distributions of vegetation exist. Shnerb and collaborators 
utilized a hybrid approach to this problem [13] . They considered the deterministic equa- 
tions governing the vegetation model, but introduced an additional "integration trick" 
aimed at incorporating in some effective way the discrete nature of the underlying individ- 
uals/plants and, hence, to incorporate indirectly stochastic demographic effects. Thus, 
the question arises as to whether it is possible to account for these discreteness-induced 
effects purely at the coarser level of the population, that is, introducing stochasticity in 
the deterministic equations in a more systematic and controlled way. Again, thanks to 
the methods and concepts developed in field theory [11], we know that the answer is 
yes. As in the example of plankton above, it is possible to derive analytically a set of 
Langevin equations starting from the discrete model by using the same approach. Once 
more, however, the resulting set of equations is technically difficult to deal with. 

A third example to be discussed is the formation of regular, almost crystal-like, pat- 
terns in bacterial colonies. To this end, we use the simple "Brownian bug model", in 
which random walkers (bugs) diffuse, branch, and die at some rates. While deterministic 
approaches to this model lead to too-perfect orderings, incompatible with fundamental 
principles of physics (see below), the inclusion of stochasticity leads to much more realistic 
patterns, as well as accurate and precise predictions [T5l [161 EE] ■ 

Summing up, at the end of this discussion, we are left with the following dichotomy: 
either we have i) an individual-based description very detailed but lacking emphasis on 
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large-scale features and, hence, handicapping the understanding of the emerging phe- 
nomenology, or ii) a deterministic global description in terms of density fields, which 
emphasizes large-scale properties but potentially misses important features, specially in 
low- dimensional systems and for low densities. As described above, in some cases a third 
way exists: it is possible, starting from individual-based models, to derive analytically 
stochastic continuous (Langevin) equations. They usually are extensions to their deter- 
ministic counterparts that include additional demographic-noise terms. For instance, a 
recent example of the successful application of this systematic approach to the spatial 
Lotka-Volterra model can be found in [18]. Unfortunately, such Langevin equations are 
usually difficult to treat through computational studies, leaving analytical approaches 
(which may or may not be feasible) as the only available options [19] . 

The essential problem one faces in trying to study computationally Langevin equations 
with demographic noise is that standard integration methods generate unrealistic negative 
density values [20]. Demographic noises are proportional to the square- root of the density 
field and, owing to this, whenever the local values of the density become close to zero, 
the noise term is much larger in magnitude than any other term in the equation. This, 
combined with the random sign of the noise term, leads standard integration schemes 
ineluctably to negative densities and, hence, to numerical instabilities. However, a novel 
and powerful integration method has been developed to deal with demographic noise in 
an accurate and precise way [2"T[ |2"2"] . The method is a split-step algorithm in which the 
system is discretized in space, and in which temporal integration is implemented in two 
steps: (i) The noise term plus linear terms in the deterministic dynamics are treated in an 
exact way by sampling the conditional probability distribution coming out of the (exactly 
solvable) associated Fokker-Planck equation. By sampling such a distribution, an output 
is produced at each site, (ii) Then, the remaining deterministic terms are integrated 
using any standard scheme, choosing as initial conditions the output of the previous step 
at each site. This scheme is able to avoid the difficulties associated with negative density 
values, and converges to the Langevin equation solution for discretization times that need 
not be as small as in the standard integration approaches (see [2T[ [22] for further details). 

In this paper, we apply the split-step integration scheme to Langevin equations de- 
scribing the three examples discussed above. We focus on the (hardest-to-analyze) small- 
density limit, i.e. close to extinction. To avoid further difficulties (stemming mostly from 
noise cross-correlations and conserved noise terms) we use a minimal description in which 
we take the deterministic equations for each case and add the simplest possible form of 
(biologically reasonable) demographic noise to each of them. This can be considered as 
either a stochastic extension of the deterministic models, or as a simplification of the more 
complex Langevin equations analytically derived from the corresponding individual-based 
model, in which higher-order irrelevant terms are omitted. By using the split-step inte- 
gration scheme, we scrutinize numerically the role of demographic noise in the generation 
of patterns in the stochastic versions of the examples presented above. Observe that, in 
all these examples, when extinction is reached (i.e. when the density of phytoplankton, 
vegetation or bacteria vanishes) the system remains trapped in such a state indefinitely: 
phytoplankton, vegetation, and bacteria do not arise spontaneously and, therefore, the 
empty state is an absorbing state. Thus, it is also interesting to investigate if the behavior 
of our Langevin equations belongs to any known universality class for systems with ab- 
sorbing states [231 [2U 1213 I2E]- We show that this minimal continuous stochastic approach 
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keeps the relevant fluctuations present at the discrete level in the studied examples and, 
in consequence, is able to reproduce the most significant phenomenology and key features 
observed in the microscopic-level counterparts. 



2 Patterns in Oceanic Plankton 

The Levin-Segel model (LSM) is defined by a set of two deterministic equations accounting 
for the basic interactions between a population of phytoplankton (autotrophic organisms) 
and zooplankton (heterotrophic organisms that feed on phytoplankton, i.e. grazers) [51 E]. 
If p(x, t) is the density of phytoplankton at position x and time t, and 0(x, t) the density 
of zooplankton at the same coordinates, the LSM may be written as: 

T QM . / 9 t p(x, t) =ap + bp 2 - wfxj) + DV 2 p , v 

1 <9 t 0(x, t) = w 2 p<P - A0 2 + D 2 V 2 <P, [i) 

where a, b, w, D, w 2 , X and D 2 are constants [27]. The first two terms on the right hand 
side of the first equation describe the growth and replication of phytoplankton cells, the 
third one the mortality due to zooplankton consumption, and the fourth one the diffusion 
of cells. Similarly, the first and last terms on the r.h.s. of the second equation describe the 
growth of the zooplankton population due to grazing and effective diffusion, respectively, 
while the second one is a saturation term imposing a certain carrying capacity. 

Together with the trivial extinction state, p(x, t) = 4>(x, t) = 0, and the (biologically 
implausible) state in which <fi vanishes and p diverges, Eqs.Q have a spatially uniform 
stable equilibrium in which p = aX/(ww 2 — bX) and = a«j 2 /(™2 — bX), provided 
ww 2 > bX and w 2 > b. In the usual way, the separation of scales between diffusion 
constants allows for Turing patterns to emerge [I]; in particular, when the ratio D 2 /D 
is larger than a certain constant K cr (which is a combination of some parameters [7J), 
the homogeneous state becomes unstable, entailing pattern formation. Thus, the line 
D 2 /D = K cr > 1 separates the zones of the parameter space where it is possible to 
generate Turing patterns from those where the homogeneous state is stable (see blue 
curve of Fig. 1 in [H] ) . 

As mentioned above, an individual-level version of that model was introduced by 
Butler and Goldenfeld [9]. Based on the interactions between phytoplankton (P)and 
zooplankton (Z) proposed by the LSM, it explicitly writes reaction equations for the 
different sources of reproduction and mortality: 

a b/V 
P > PP PP PPP . . 

w/V k/V v I 

PH HE HH -U H 

(where V is the volume of a well-mixed patch where these reactions are considered, and 
w = w 2 is assumed), together with the passive diffusion of the agents. From the proposed 
set of individual interactions, by using standard field-theory techniques [TTj EH] , a set of 
coupled Langevin equations for the density fields was deduced. Such equations coincide at 
the deterministic level with Eq.flTJ, but they include additional noise terms (see Eq.(15) 
in [9]). As shown analytically by the authors, this stochasticity greatly enlarges -owing to 
a resonant amplification of fluctuations [29]- the region in which patterns emerge, leaving 
only a small region of parameter space with homogeneous stable solutions. 
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The actual complex form of the noise in these Langevin equations makes it difficult 
(if not impossible) to verify numerically the analytical result at the mesoscopic level. Of 
course, it is always possible to perform direct simulations at the microscopic level, but 
here we are interested in the mesoscopic -Langevin- description. In order to incorporate 
to the relevant individual-level fluctuations to the set of discrete equations in Eqs.flTJ), we 
add demographic-noise terms to each of the differential equations in the simplest possible 
form compatible with biological constraints: 

o t qi\/r . f <V(X, t) =ap + bp 2 - Wp(j) + DV 2 p + (Ty/p 7?(x, t) , . 

1 <9 t 0(x, t) = w 2 p(t> - A0 2 + D 2 V 2 + a 2V /p0 %(x, t), {6> 



where a and a 2 are constants, and rj and 7] 2 are delta-correlated Gaussian noise (i.e. 
white noise) terms [30J. Heuristically, one can argue that -as it is usually the case in 
particle systems, and as a direct consequence of the central limit theorem- the noise 
amplitude has to be proportional to the square-root of the involved densities. Observe 
that, in the case of zooplankton, reproduction and mortality are related to the presence 
of phytoplankton; therefore, the variance of the new noise term must be proportional 
to both zooplankton and phytoplankton densities; for phytoplankton, the leading order 
of approximation leaves the variance of the noise proportional to the phytoplankton 
density. Alternatively, it is not difficult to show that these terms, present in the detailed 
derivation in [9], are the leading terms in that approach. Higher-order corrections are 
irrelevant in the renormalization group sense and noise cross-correlations, which might not 
be irrelevant, are anyhow perturbatively generated. The advantage of our simplification 
is that, now, the resulting set of Langevin equations can be numerically studied by direct 
application of the split-step integration scheme for Langevin equations with square-root 
noise [2T|l22]. 

We use the integration scheme in a 2-dimensional lattice of lateral size L. We study 
the change in time of the spatially-averaged density of both fields, p and 0, starting from 
a homogeneous state, and determine how that behavior is altered by changing the control 
parameter a while fixing the rest of parameters in Eqs.(!51). As happens with the original 
LSM, we find two phases: an active phase, where both phytoplankton and zooplankton 
reach a stationary non-trivial state, and an absorbing phase, where the primary extinction 
of phytoplankton leads unavoidably to zooplankton extinction. The integration algorithm 
permits us to reproduce generically patterns almost all along the active phase, allowing 
for a numerical verification, at the mesoscopic level, of the analytical findings in [9j| (see 
FigHD. 

Note that the lack of a saturation term for phytoplankton in Eqs.([3]) could lead to 
the erroneous conclusion that its averaged density could grow boundlessly in the su- 
percritical phase (i.e. no stationary state would be possible). Interestingly, it is the 
coupling between fields in the zooplankton equation that keeps the phytoplankton field 
bounded: when phytoplankton density grows, the positive coupling in the equation 
entails a consequent zooplankton density growth which, in turn, has a bounding effect 
on phytoplankton density. Thus, zooplankton acts as an inhibitor and stabilizes the 
phytoplankton population. 
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Figure 1: Up: Generic phase diagram of the stochastic version of the LSM. Down: Representative 
examples of the disordered patterns formed by the phytoplankton population in the supercritical phase. 

Critical properties and universality 

To study universality issues at the transition point, and in deference to the simplest 
biology (without active aggregation), we relax the assumption that D 2 S> D and, actually, 
set D = D 2 to a small value -see caption of FigfSJ This enhances the convergence to 
the asymptotic state while avoiding the presence of "Turing-like" effects (as D/D 2 = 1). 
As the condition b ^ in Eqs.([T]) is essential for both the original [7] and extended [9] 
versions of the model to show a non-homogeneous solution, we also impose b = in order 
to challenge the ability of our new description to show patterns. 

As can be seen in the left panel of FigfS] (inset), starting from a homogeneous initial 
condition, for values of a below a c = 0.31740(5) and large system sizes there is an initial 
power-law decay of the phytoplankton population density which eventually goes extinct 
in a finite time. For values of a above such value, the population eventually reaches a 
stationary state with a non-zero density for both species. At a = a c , the system undergoes 
a second-order phase transition characterized by the scale-invariant decay of the density 
following a power law of exponent 9 = 0.49(5). 

As a consequence of the existence of a critical point, one expects scale-invariant be- 
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Figure 2: Dynamic behavior of the S-LSM with parameters b = 0, w = W2 = 1, D = D2 = 0.25, 
A = 1, a = 02 = \/2 and dt = 0.1. Left: Temporal decay (in MonteCarlo steps) of the spatial average 
of phytoplankton density, starting from a homogeneous distribution, for different system sizes. Inset: 
Resulting curve of multiplying ~p(t) by t e , using a system of linear size L = 2 11 (i.e. 4, 194, 304 sites); the 
horizontal curve corresponds to the critical point. Right: Mean quadratic radius (upper curve), number 
of sites with non-zero density and survival probability (lower curve) for the spreading of a localized seed 
of phytoplankton and zooplankton in an otherwise empty system. 



havior of other quantities, i) Performing simulations starting from an initial small density 
for both fields, localized at a single point, we can study how the population spreads over 
an otherwise empty system. At the critical point, the mean quadratic radius of the pop- 
ulation (mean square distance of the population border to the original "seed"), r 2 , the 
number of lattice sites with a non-zero density of any of the two species, N s , and the 
survival probability of the population, P s , change with time following a power law of 
exponents z spr = 1.13(2), r\ = 0.23(2) and 5 = 0.42(2), respectively (see right panel of 
FigJSj). ii) The stationary density of the smaller sizes (see left panel of FigJ2J) on the linear 
system size L, follows a power law described by p st ~ L 13 ^^. This is also the case of the 
survival time of experiments started from homogeneous conditions, t sur , which depends 
on L following a power law of exponent v\\/v± (see table 1 for actual values, and left panel 
ofFigE]). 
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0.4505(10) 


0.4505(10) 


0.2295(10) 


1.1325(10) 


0.795(4) 


1.766(2) 


S-LSM 


0.49(5) 


0.42(2) 


0.23(2) 


1.13(2) 


0.83(5) 


1.73(5) 


S-Shnerb 


0.47(5) 


0.48(2) 


0.23(2) 


1.16(2) 


0.79(5) 


1.75(5) 



Table 1: Critical exponents obtained with the simple stochastic versions of the determinis- 
tic models analyzed in the text, plus the corresponding values for the directed percolation 
(DP) universality class. The exponents for the "Brownian Bug" model, in agreement with 
the DP universality class, can be found in [45] . 

The right panel of Figj3] is an illustrative example of the spatial distribution of the 
phytoplankton population in the supercritical stationary states. As we can see, phyto- 
plankton form disordered patterns with patches of different sizes with non-zero density. 

In summary, the set of equations Eqs.flH]) with 6 = and D / D 2 = 1 is not only able to 
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Figure 3: Up: Finite-size scaling of the stationary values depicted in Fig[2](left panel) and the survival 
time for homogeneous experiments when using a threshold for the survival probability of 0.1. Down: 
Typical snapshot of the phytoplankton density in the supercritical phase. 



describe patchiness in its spatial distribution (a common feature or real phytoplankton- 
zooplankton populations [7J EI]), but also shows a richer, non-trivial phenomenology not 
observed in the original deterministic formulation, which stems from stochasticity. Care- 
ful inspection reveals that the derived set of critical exponents (see table 1) corresponds 
to the universality class of directed percolation, which controls many transitions into ab- 
sorbing states [23"j EH EI3 EE] ■ The reasons why this model, as well as the rest of models 
discussed in the paper, belongs to the Directed Percolation universality class are dis- 
cussed in the Appendix. Analogously, our approach also permits us to study a region of 
parameters in which either supercritical DP clusters or other non-trivial patterns emerge 
(see FigJT]), and verify that this region is indeed largely enhanced with respect to the 
deterministic predictions. 



3 Patterns in Vegetation in Semi- Arid Environments 

Aimed at describing the patterns of vegetation observed in semi-arid zones, Klausmeier 
introduced a model that keeps track of the interactions between vegetation biomass and 
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water in such ecosystems [12J. This model, as well as an extension to it including a 
more complex water diffusion term [13], leads to pattern formation only in the "Turing 
limit" (i.e. very different diffusion rates for the vegetation and water density fields), or 
in inhomogeneous media. 

Shortly afterward, Shnerb et al. revisited the problem in a series of papers [HI [321 E31 
|3l], introducing a simplified, combined version of the two previously mentioned models 
with the intention of reproducing the disordered patterns frequently observed in semi-arid 
plains. If p represents the biomass density, and water density, the equations read: 



where a, R, w, u>2, k, and D are constant parameters, and v is a constant velocity 
vector. In this model, vegetation (first equation) grows in the presence of water, there 
is an effective diffusion due to, e.g. seed dispersal, and there is also competition for 
water, here represented by the first term on the r.h.s. On the other hand, water density 
(second equation) increases due to precipitation (R) and decreases due to consumption 
and evaporation (second and third term on the r.h.s.); it can also flow from high to low 
places, as represented by the last term of the r.h.s. 

Eqs.fll]) exhibit two different homogeneous equilibria: a trivial state where vegetation 
goes extinct and water reaches the stationary value imposed by precipitation (p = 0, <fi = 
R), and a nontrivial one where both reach a stationary state given by (p = (R — a)/(aw 2 ), 
cf) = a). In the absence of cross- diffusion terms like the one proposed in [T5] or anisotropies 
such as the one imposed by v ^ 0, the system lacks unstable homogeneous solutions that 
can be perturbed in order to obtain the desired patterns. 

Shnerb and coauthors realized that Eqs.fll]) are able to show realistic disordered pat- 
terns when a seasonal removal of biomass below a certain threshold is introduced [H"| 132] , 
and/or an additional selective mortality that depends on the neighboring biomass [T3] . 
Both mechanisms account in some effective way for stochastic effects at low densities, 
although they somehow mix different scales of description. 

We follow now the steps of the previous section to reproduce disordered patterns in 
a systematic way without resorting to integration prescriptions [35] . In analogy with the 
previous section, we introduce a simple demographic noise term in the equation for the 
vegetation biomass density. As in the case above, we set v = (i.e. no anisotropy or 
inhomogeneous soil) and D/D 2 = 1. The resulting equations for the stochastic version 
of the Shnerb model (S-Shnerb) are: 



where the terms and constants are as explained for Eq. p]) , and the new term is identical 
to the one introduced for phytoplankton in the previous example. 

Using again the split-step integration scheme, and choosing the precipitation density 
R as tuning parameter, we observe that, for R < R c = 0.40835(5), the spatially-averaged 
density of vegetation biomass decays until the population becomes extinct, while for 
R > R c , it reaches a non-trivial stationary state. The system undergoes a phase transition 
between these two states at R = R c , where p decays asymptotically toward extinction 
following a power-law. On the other hand, the averaged water field, 0, eventually reaches 



Shnerb 




d t p(yi,t) = —ap + wptfi + DV 2 p 

<9 t 0(x, t) = R - w 2 p(p - + D 2 V 2 (p - v. V0, 



(4) 



S- Shnerb 
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a stationary value in any of the three situations. For R < R c , <p grows with R, while for 
R > R c , that value decreases as R (and p st ) increases. 

Once more, the biological processes involved shed light on how the stationary state 
is achieved in the absence of a saturation term in Eq.(jS]). As described above, when the 
vegetation population goes extinct, larger values of the precipitation constant R lead to 
larger values of the average water density 0. On the other hand, there is a change in 
trend when vegetation survives. For values of R into the supercritical phase, the larger 
the precipitation parameter is, the better the conditions for vegetation growth are. In 
consequence, the average biomass p reaches a larger value at the stationary state, and the 
improved consumption influences the average value of the water density, which decreases 
with the increase of R. Therefore, increasing R entails a decrease in the average of linear 
terms in the p equation, which keeps bounded the biomass field. The competition for 
the available resource (water) acts as an inhibitor, which allows the vegetation density to 
reach a well-defined stationary state, acting as an effective carrying-capacity term. 
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Figure 4: Behavior of Eqs.© with a = 0.2, w = 1, w 2 = 1.2, re = 1, D = D 2 = 0.25 and a = s/2. Up: 
Spreading observables (main panel) and stationary observables (inset) for R = R c = 0.40835. Down: 
Typical snapshot of the vegetation density in the supercritical phase. 

Analyzing both the dynamic and stationary behaviors of Eqs.flH]), we obtain the curves 
and exponents shown in Fig H] and table 1. As in the previous case, the existence of a 
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critical point entails scale invariance that is translated into power-law observables and 
disordered clusters of vegetation of any size. For R = R c and above, it is possible to find 
the desired disordered patterns with these equations (see right panel of FigHJ). As in the 
example above, the measured critical exponents agree (within error bars) with those of 
the DP class (see table 1 and Appendix). 



4 Brownian bugs 

As a third and final example, we present an individual-based model: the "interacting 
Brownian bug" (IBB) model [To! [TB] . It consists of branching- annihilating Brownian 
particles (bugs, bacteria, etc) interacting with each other within a finite distance, / [T5] . 
Particles move off-lattice in a (^-dimensional space, and their dynamics is such that they 
can diffuse at rate 1, performing Gaussian random jumps of variance 2D; disappear 
spontaneously, at some rate (3q; or branch, creating an offspring at their location with a 
density-dependent rate A modeling competition: 

A(j) = max{0, A - Nt(j)/N sat }, (6) 

where j is the particle label, Ao (reproduction rate in isolation) and N sat (saturation 
number) are fixed parameters, and Ni(j) stands for the number of particles within a 
radius / from j. The control parameter is given by \l = Ao — A> the system is active 
above some \i c and absorbing below. In the active phase, owing to the density- dependent 
dynamical rules, particles group together forming clusters with a characteristic size. Well 
inside the active phase, when these clusters start filling the available space, they self- 
organize in spatial patterns with remarkable hexagonal order. 

The IBB model can be cast into a continuous stochastic equation [TSJ [TT] . Indeed, by 
applying the same standard techniques as above, a Langevin equation for the density p 
of bugs can be derived [UJ : 

^^-=pp + DV 2 p-^- [ dyp(y,t) + a^(x,t), (7) 

01 iV S J\ X -y\<R 

where the noise amplitude a is a function of the microscopic parameters, r/(x, i) is a 
normalized Gaussian white noise, and higher order terms and cross-correlations (irrelevant 
in the renormalization group sense) have been neglected. Remarkably, the deterministic 
part of Eq.flTJ, including a non-local saturation term, is a particular case of the equation 
proposed by Fuentes et al. |36j as a model for competition-induced pattern formation, in 
which the non-linear term is: 



p(x,t) j dyF{x-y)p{y,t) 



and F is a generic kernel or influence function (which becomes a step function in Eq.flTJ). 

A numerical integration of Eq.flTJ relying on the split-step scheme has been recently 
performed [UJ. It reproduces very accurately the phenomenology of the microscopic 
IBB model described above: there is a critical point in the DP class separating the 
patterned active phase from the absorbing phase. Two important observations are in 
order: i) As carefully discussed in [TT] and illustrated in FigJS] the resulting ordering 
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Figure 5: Different patterns for different parameter values as obtained from a numerical integration of 
Eq. ([7]). The amplitude of the noise decreases monotonously from (a) to (d). 

is not "perfect" as would be the case for the corresponding deterministic equations. 
The presence of demographic noise induces "defects" and the resulting distribution of 
clusters is not a perfect regular crystal-like structure (which is actually precluded in two 
dimensional equilibrium systems by fundamental physics principles, i.e. the Mermin- 
Wagner theorem). In fact, as shown in [T7| the spatial structure (including defects, as 
well as their dynamics and statistics) is perfectly described by the theory of (equilibrium) 
two-dimensional melting, ii) As shown in [37] the case in which the kernel is Gaussian 
does not lead generically to clustering but rather to homogeneous solutions. However, 
it has been recently observed that, by introducing demographic noise, robust patterns 
similar to those reported above emerge. These are purely noise-induced patterns [3"5| I3"5]. 
that can be reproduced by our approach. 

5 Conclusions 

Connecting scales of observation is an important and delicate task. Important features 
of the individual level can be lost during the scaling up to the level of the population 
unless a rigorous approach is followed. In cases where a careful deduction of coarse- 
grained equations can be performed, these can be difficult to deal with. A trade-off 
between realism and complexity is to be resolved attending to the desired features to be 
reproduced by the continuous equations. 

Simple forms of demographic noise, added to standard deterministic equations for 
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different problems involving pattern formation, suffice to provide a much more precise 
description of the underlying "microscopic" dynamics. We have considered three different 
examples: plankton dynamics, a model for vegetation growth in semiarid environments, 
and a model for interacting and diffusing particles. In all these cases, by adding de- 
mographic noise -even in its simplest form- the regions in the parameter space where 
patterns appear are greatly enlarged. Intrinsic noise also introduces a non-trivial phase 
transition separating the survival state from the absorbing or extinction one in all the 
reported examples. 

This conclusion is backed by computational studies of the resulting stochastic/Langevin 
equation for each of the examples. Numerical integration of Langevin equations with 
demographic ( "square-root" ) noise is feasible owing to a recently proposed split-step in- 
tegration scheme, avoiding the otherwise unavoidable integration instability associated 
with this kind of noise. 

By using this split-step scheme, we have observed patchy distributions in regions for 
which the corresponding deterministic approach would only lead to homogeneous steady 
states. We have located a non-trivial phase transition separating the patchy active phase 
from the extinction/absorbing phase. Moreover, we have measured the critical exponents 
associated with such phase transition, finding that the universality class of the investi- 
gated examples is directed percolation. This universality class is paradigmatic of systems 
with absorbing states in the absence of other relevant symmetries. In the two first dis- 
cussed examples, we have seen that the coupling of activity to a non-conserved, diffusing 
field does not alter the universality class. Moreover, the absence of an explicit saturation 
term for activity results to be irrelevant here, as the (inhibitory) ecological interactions 
between the fields suffice to maintain the density of the interacting species bounded into 
finite intervals, allowing for a well-defined stationary state (either coexistence or extinc- 
tion). 

In summary, properly derived Langevin equations -constructed either as an extension 
of deterministic models with additional demographic noise terms or as simplifications 
of more formally derived complicated stochastic equations- provide a highly valuable 
tool for the analysis of population-level ecological properties. This is particularly true 
as long as numerical studies of such Langevin equations (a direct way to explore their 
phenomenology) are feasible. Simple Langevin equations combined with the split-step 
integration scheme provide us with a powerful tool to scrutinize population-level features 
keeping the relevant effects stemming from the underlying discrete nature of individuals. 
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Appendix 



We have numerically shown that the critical exponents measured in all of the exam- 
ples discussed in this paper are directed-percolation like. Here we justify that result by 
analyzing their corresponding Langevin equations Eqs.flH]), Eqs.^SJ), and Eq. ([7]). 

As conjectured years ago by Janssen and Grassberger [ID] all systems exhibiting a 
phase transition into a unique absorbing state, with a single-component order parameter, 
and no extra symmetry or conservation law, belong to the same universality class whose 
most distinguished representative is directed percolation (DP). In a field theoretical de- 
scription this universality class is represented by the Reggeon Field theory (RFT), [10] 
which in terms of a Langevin equation reads: 



where p(x, t) is the density field at position x and time t, a and b are parameters, and 
t](x,t) is a delta-correlated Gaussian noise, < f](x,t)i](x' ,t') >= D5(x — x')5(t — t'). 

The previous conjecture was confirmed in a large number of computer simulations, 
series expansion analysis, field theoretical studies, etc. Indeed, the DP universality class 
has proved to be extremely robust against the modification of many details in the mi- 
croscopic models. The conjecture of universality has been extended systems with an 
infinite number of absorbing states |41| . Also, Grinstein et al. extended the conjecture 
to the case of multicomponent systems [12]. Considering, for the sake of simplicity a 
two-component system with absorbing states, it can generically be described by a set 
of Langevin equations whose linearized dynamics (ignoring the Laplacian terms and the 
noise) is generically: 

d t pi(x,t) = a 1A px + a h2 p2 

d t p 2 (x,t) = a 2 ,ipi + a 2j2 p 2 ' ' 

and diagonalizing the matrix of the linear coefficients: 



where Ai and A 2 are the associated eigenvalues and £i and £ 2 are the corresponding 
eigenvectors. The existence of the absorbing state implies that the two eigenvalues are 
negative. At the critical point one of the two eigenvalues, say Ai vanishes. If the other one 
does not, it remains negative, and then fluctuations in £ 2 continue to decay exponentially 
with time. This implies that £ 2 does not experience critical fluctuations in the vicinity of 
the transition, and so can be integrated out of the problem. Hence, asymptotically, the 
set of Langevin equations can be reduced to single relevant equation which is precisely 
Eq.®, ensuing DP behavior. 

On the other hand, if the secondary field is conserved (which, in particular, enforces 
a 2j i and a 2 2 to vanish) or there is a hidden symmetry imposing Ai and A 2 to vanish at 
the same point, this conclusion can break down. In both of these cases, at the system 
critical point, the two fields are simultaneously critical (or "massless" or "gapless" using 
the field theory jargon). 




(9) 




(11) 
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Systems in which the density field is coupled to a secondary conserved field indeed 
exhibit non-DP critical behavior. This is the case of, for instance, i) systems related 
to self-organized criticality in which the activity field is coupled to a conserved non- 
diffusive background field [13] and ii) systems in which the secondary field is conserved 
but diffusive [lUSS]. 

Let us now return to the three models discussed in this paper. 

a) Let us first scrutinize the simpler case of vegetation patterns as described by Eqs.([5]). 
At the critical point the expectation value of p vanishes, while converges on 
average to a value Rj k and, for generic values of p, it converges to st — (R — 
W2 p)/k- Let us remark that these are the noiseless or mean-field expectation 
values. Defining a new field ip = — <p st and linearizing the dynamics as above, it is 
straightforward to see that this is an instance of the non-symmetric case discussed 
above in which the linear terms in both equations do not vanish simultaneously: 
DP behavior is predicted. Observe also that, after shifting the secondary field, a 
standard saturation term proportional to — p 2 appears in the first equation, i.e. the 
competition for the available resource (water) acts as an inhibitor generating the 
necessary effective carrying capacity term (see main text). 

b) The case of the stochastic Levin-Segel model, as defined by Eqs.tjl]) is a bit trickier. 
Observe that on average one expects (at a noiseless or mean-field level) = pw 2 /X, 
which implies that when p vanishes so does opening the door to non-DP scaling. 
However, closer inspection reveals that the noise in the equation for is a higher or- 
der one as it is proportional to \f p 0. It is therefore irrelevant in the renormalization 
group sense around the DP-like fixed point (indeed, we have verified numerically 
that the critical exponents of the modified version of Eqs. (JH]) without such a term 
coincide, within error bars, with those of the original ones -results not shown). In 
this way, the equation for dt<j) becomes asymptotically a deterministic one. Then, 
0(x) can be written as a function of p(x,t), i.e. 4>(x,t) = F(p(x,t)) where F is an 
unspecified function which can be formally obtained by integration of the equation 
for <9 t 0. Expanding such a function in power series of p, plugging in the result in the 
equation for dtp and neglecting higher-order terms, one readily recovers Eq. (Q, 
and hence DP scaling. The field is just a slave mode of p (proportional to it up 
to leading order), inheriting all its critical properties from it. 

c) Finally, for the non-local Langevin model for Brownian bugs, it is easy to see that 
by performing a coarse graining in which scales below / are scaled-up, one recovers 
effectively Eq.(|9]). The main point is that nonlocal but finite-range interaction 
becomes local in the renormalization group perspective once a sufficient amount of 
coarse graining has been considered. See [T7] for further details. 
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